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Abstract. In this paper we present two optimized eight-step symmetric implicit 
methods with phase- lag order ten and infinite (phase-fitted) . The methods are con- 
structed to solve numerically the radial time-independent Schrodinger equation with 
the use of the Woods-Saxon potential. They can also be used to integrate related 
IVPs with oscillating solutions such as orbital problems. We compare the two new 
methods to some recently constructed optimized methods from the literature. We 
measure the efficiency of the methods and conclude that the new method with 
infinite order of phase-lag is the most efficient of all the compared methods and for 
all the problems solved. 

Keywords: Schrodinger equation, orbital problems, phase-lag, initial value prob- 
lems, oscillating solution, symmetric, multistep, implicit 

PACS: 0.260, 95.10.E 



1. Introduction 

The radial Schrodinger equation can be written as: 



y"{r)={-^-j^ + V{r)-E]y{r) (1) 



r 



where ^^^^ is the centrifugal potential, V{r) is the potential, E is the 

energy and W{r) = ^^^3^ + V{r) is the effective potential. It is vahd 
that hm Vir) = and therefore hm Wir) = 0. 

r— >oo r— >oo 

We consider E > and divide [0, 00) into subintervals [oj, bi] so that 
W{r) is a constant with value Wi. After this the problem (1) can be 
expressed by the approximation: 

y'^ = iW-E)y„ (2) 
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whose solution is: 

yi{r) = Ai exp {\/w - E r) + Bi exp (- ^/w - E r) , 
A,, B, G R. 

Many numerical methods have been developed for the efficient so- 
lution of the Schrodinger equation and related problems. For exam- 
ple Raptis and Allison have developed a two-step exponentially-fitted 
method of order four in [5]. More recently Kalogiratou and Simos have 
constructed a two-step P-stable exponentially-fitted method of order 
four in [6]. 

Some other notable multistep methods for the numerical solution 
of oscillating IVPs have been developed by Chawla and Rao in [3], 
who produced a three-stage, two-Step P-stable method with minimal 
phase-lag and order six and by Henrici in [4] , who produced a four-step 
symmetric method of order six. Also Anastassi and Simos have devel- 
oped trigonometrically fitted six-step symmetric methods in [108] and a 
six-step P-stable trigonomctrically-fittcd method in [110] and Panopou- 
los, Anastassi and Simos have constructed two optimized eight-step 
symmetric methods. 

Also some research work in numerical methods can be found in [1]- 
[137]. 



2. Phase-lag analysis of symmetric multistep methods 

For the numerical solution of the initial value problem 

y" = f{x,y) (4) 

multistep methods of the form 

m m 

X! ^iVn+i = h^J2 ^ifi^n+i, Vn+i) (5) 
i=0 j=0 

with m steps can be used over the equally spaced intervals {a^jj^Q G 
[a, h] and h = \xi^i — Xi\, i = 0(l)m — 1. 

If the method is symmetric then Oj = am-i and bi = bm-i, i = 

o(i)LfJ- 

Method (5) is associated with the operator 

m m 

L{x) = ^ aiu{x + ih) -h'^^ biu"{x + ih) (6) 

i=0 i=0 

where n G C^. 
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DEFINITION 1. The multistep method (5) is called algebraic of order 
p if the associated linear operator L vanishes for any linear combination 
of the linearly independent functions 1, 

When a symmetric 2A;-step method, that is for i = — is applied 

to the scalar test equation 

y" = -^'y (7) 

a difference equation of the form 

Ak{v)yn+k + ••• + Ai{v)yn+i + Ao{v)yn 

+Ai{v)yn^i + ... + Ak{v)yn-k = (8) 

is obtained, where v = coh, h is the step length and Aq(v), Ai(v), . . ., 
Ak{v) are polynomials of v. 

The characteristic equation associated with (8) is 

Ak{v)s'' + ... + Ai{v)s + Ao{v) + Ai{v)s~'^ + ... + Ak{v)s-'' = (9) 



Prom Lambert and Watson (1976) we have the following definitions: 

DEFINITION 2. A symmetric 2fc-step method with characteristic equa- 
tion given by (9) is said to have an interval of periodicity (0, Vq ) if, for 
all V G {0,Vq), the roots Si,i = 1(1)2A; of Eq. (9) satisfy: 

si = e^^(^), S2 = e-^^(^\ and \si\ < 1, i = 3(1)2A; (10) 

where 9{y) is a real function of v. 

DEFINITION 3. For any method corresponding to the characteris- 
tic equation (9) the phase-lag is defined as the leading term in the 
expansion of 

t = v-e{v) (11) 

Then if the quantity t = 0{v'^~^^) as v ^ oo, the order phase-lag is q. 

THEOREM 1. [1] The symmetric 2k-step method with characteristic 
equation given by (9) has phase-lag order q and phase-lag constant c 
given by 



-cvi+^+0{vi+^) = 



2Ak{v) cos{kv) -\- ... -I- 2Aj{v) cos{jv) -\- ... -\- Ao{v) 
2k'^Ak{v) + ... + 2fAj{v) + ... + 2Ai{v) 

(12) 
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The formula proposed from the above theorem gives us a direct 
method to calculate the phase-lag of any symmetric 2k- step method. 

In our case, the symmetric 8-step method has phase-lag order q and 
phase-lag constant c given by: 

-cv'}+^ + 0{v'i+^) = 

2A4{v) cos(4i;)+2A3 {v) cos(3t))+2A2 jv) cos(2d)+2Ai (v) cos{v)+Ao (v) (13) 
32A4(?;)+18A3 (v)+8A2 {v)+2Ai (v) 



3. Construction of the new optimized multistep methods 

We consider the eight-step symmetric implicit methods of the form: 

Vi = -y-A - 03(2/3 + y-3) - 02(2/2 + y-2) - ai{yi + y-i) 

(64(/4 + /-4) + 63(73 + /-3) + 62(/2 + f-2) + 61 (/l + /-l) + 60/0) 

(14) 

where 03 = —2, 02 = 2, ai = —1, 

yi = y{x + ih) and /j = f{x + ih, y{x + ih)) 



3.1. First optimized method with infinite order of 

PHASE-LAG (phase-fitted) 

Wc want the first method to have infinite order of phase-lag, that is 
the phase-lag will be nullified using 64 coefficient. 

We satisfy as many algebraic equations as possible, but we keep 64 
free. After achieving 10th algebraic order, the coefficients now depend 
on 64: 



(15) 



60 = 7064-1^, 61 = -5664 + 1^, 

62 = 28 64 - |§f|, 63 = -8 64 + 
and the phase-lag becomes: 

-P^ = TWS' where 

A = 24192 (cos {v)f + 24192 (cos {v)f v%4^ + 17671 (cos {v)f v'^ 
-96768 (cos (v)f v^b^ - 24192 (cos {v)f + 14152 (cos {v)f v%i 
-12096 (cos {v)f - 11811 (cos {v)f + 2109 cos {v) v^^ 
15120 cos [v) - 96768 cos {v) v%4. - 409 v'^ + 24192 v'^b4^ - 3024 and 
B = 12 + 25v'^ 

(16) 
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SO by satisfying PL = 0, we derive 



1 c 



where 



"4 — 24192 D ' 

C = 24192 (cos {v)f + (17671 v'^ - 24192) (cos {v)f 

- (12096 + 11811t;2) (cos (w))^+ (17) 

(15120 + 2109^2) cos {v) - 409 i;^ _ 3024 

D = v^ {cos iv)^ - 4 cos {vf + 6cos {vf - 4 cos {v) + l) 
where v = u h, u is the frequency and h is the step length used. 

3.2. Second optimized method with tenth order of 

phase-lag 

For this method we use all bi coefficients for achieving maxiniTim al- 
gebraic order or maximum phase-lag order. After achieving maximum 
algebraic order, that is ten, the coefficients become: 

_ 17273 _ 280997 33961 _ 173531 45767 

° ~ 72576' ^ ~ 181440' ^ ~ ~ 181440' ^ ~ 181440' ^ ~ 725760 ' 



If we repeat the procedure of the previous section and expand phase- 
lag using the Taylor series, we can nullify the leading term (that is 
the coefficient of h}^). However we obtain the same method as (18). 
The same method will be produced if we attempt any combination of 
algebraic order and phase-lag order. This happens due to the symmetry 
of the specific Oj. 



4. Numerical results 

4.1. The problems 

The efficiency of the two newly constructed methods will be measured 
through the integration of five initial value problems with oscillating 
solution. 

4.1.1. Orbital Problem by Franco and Palacios 

The " almost" periodic orbital problem studied by [10] can be described 

by 

y" + y = ee^^-, y(0) = l, y'{0)=i, yeC, (19) 
or equivalently by 



P_2ma.tex; 15/11/2008; 9:47; p. 5 



6 



G. A. Panopoulos, Z. A. Anastassi, T. E. Simos 



u" + u = e cos{il)x), u(0) = 1, u'(0) = 0, 
v" + v = e sm^ipx), v{0) = 0, v'{0) = l, 

where e = 0.001 and ijj = 0.01. 

The theoretical solution of the problem (19) is given below: 



(20) 



y{x) = u{x) + i v{x), u,v £ TZ 

= -^^l^^f cos(a:) + cos(^z) (21) 
= ^'i-^i' sm{x) + sin(V'x) 

The system of equations (20) has been solved for x G [0, IOOOtt]. 
The estimated frequency is = 1. 

4.1.2. Inhomogeneous Equation 

y" = -ioOy + 99 sin(t), with y{0) = 1,2/' (0) = 11, t G [0, IOOOtt]. 

Theoretical solution: y(t) = sin(i) + sin(lOi) + cos(lOt). 
Estimated frequency: w = 10. 

4.1.3. Two-Body Problem 

0, z'{0) = l,t e [0, IOOOtt]. Theoretical solution: y{t) = cos(t) and 
z{t) = sin(t). We used the estimation w = — 3- as frequency of the 

problem. 

4.1.4. Duffing Equation 

y" = -y-?/3+o.002cos(1.01t), with y(0) = 0.200426728067, y'(0) = 
0, t e [0, IOOOtt]. 

Theoretical solution: y{t) = 0.200179477536 cos(l.Olt) + 2.46946143 • 
lO-'' cos(3.03t) + 3.04014 • 10"'^ cos(5.05i) + 3.74 • 10-^° cos(7.07i) + 

Estimated frequency: w = 1. 

4.1.5. The inverse resonance problem 

We will integrate problem (1) (where r = x) with I = at the interval 
[0, 15] using the well known Woods-Saxon potential 

uq = —50, a = 0.6, xo = 7 and ui = 

a 

and with boundary condition y{0) = 0. 
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The potential V{x) decays more quickly than ^-^^7^, so for large x 
(asymptotic region) the Schrodinger equation (1) becomes 

y"{x) = (^-^^-E^y{x) (23) 

The last equation has two linearly independent solutions k x ji {k x) and 

k X ni{k x), where ji and ni are the spherical Bessel and Neumann 
functions. When x — > 00 the solution takes the asymptotic form 

y{x) Akxji{kx) — B kxni{kx) , . 

f=i D[sin{k X - TT I /2) +tan{5i) cos {kx-TT 1/2)], ^ ^ 

where 61 is called scattering phase shift and it is given by the following 
expression: 

tan (Si) = ?^(^») SiXi+i) - yjxj+i) Sjxj) 
y{xi+i) C{xi) - y{xi) C(xi+i) ' 

where S{x) = kxji{kx), C{x) = kxni{kx) and Xj < Xj+i and both 
belong to the asymptotic region. Given the energy we approximate the 
phase shift, the accurate value of which is 7r/2 for the above problem. 

We will use three different values for the energy: i) 989.701916 and 
ii) 341.495874 and iii) 163.215341. As for the frequency lo we will use 
the suggestion of Ixaru and Rizea [7] : 

/ V-E-50 X G [0, 6.5] .^.^ 
" = I xe [6.5, 15] 



4.2. The methods 

We have used several multistep methods for the integration of the 
Schrodinger equation. These are: 

— The new method with infinite order of phase-lag shown in (17) 

— The new method with eighth order of phase-lag shown in (18) 

— The P-stable method of Henrici with minimal phase-lag and order 

six [4] 

— The three-stage method of Chawla and Rao of order six [3] 

— The Classical method of Numerov 

— The P-stable exponentially-fitted method of Kalogiratou and Simos 
of order four [6] 
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— The three-step method of Adams-Moulton 



4.3. Comparison 

We present the accuracy of the tested methods expressed by the 
— log]^o(™^- error over interval) or —log^g (error at the end point), 
depending on whether we know the theoretical solution or not, ver- 
sus the log]^o(st6ps ^ stages). In Figure 1 we see the results for the 
Franco-Palacios almost periodic problem, in Figure 2 the results for 
the Inhomogeneous equation in Figure 3 the results for the Two-body 
problem and in Figure 4 the results for the Duffing equation. In Figures 
5, 6 and 7 we see the results for the Schrodinger equation for energies 
E = 163.215341, E = 341.495874 and E = 989.701916 respectively 

Among all the methods used , the new optimized method with infi- 
nite order of phase-lag was the most efficient, with the exception of the 
Duffing Equation, had almost identical results with the new method 
with phase-lag order ten. 

The difference from the other methods was about 1.2 decimal digits 
better for the Schrodinger equation for energy E = 989.701916 and 
about 0.7 d.d. for E = 163.215341 and E = 341.495874. For the other 
three problems the difference was enormous, where there was an almost 
vertical increase in the accuracy compared to the other methods. There 
were no case where the efficiency dropped below the efficiency of the 
others. 

As regards the other methods, the one of Henrici was the most 
efficient, with next the method of Chawla, the method of Numerov 
and finally the methods of Kalogiratou-Simos and Adams-Moulton. 



5. Conclusions 

We have constructed two optimized eight-step symmetric implicit meth- 
ods. The first one has phase-lag of order infinite (phase- fitted). The 
second one has phase-lag of order ten. We have applied the new methods 
along with a group of recently developed methods from the literature to 
the Schrodinger equation and related problems. We concluded that the 
new methods are highly efficient compared to other optimized methods 
which also reveals the importance of phase-lag when solving ordinary 
differential equations with oscillating solutions. 
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Figure 1. Efficiency for the Orbital Problem by Franco and Palacios 
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Figure 2. Efficiency for the Inhomogeneous equation 
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